### Re-examining Civil War: Ethnicity, Horizontal Inequality, and Space
### Grant Gordon & Evann Smith

load("imp1.RData")
load("W.RData")

ll.Wlogit <- function(par,y,W,X) {
	beta <- par[1:ncol(X)]
	rho <- par[(ncol(X) + 1)]
	inside <- 1 + exp((1 - 2*y)*(rho*W%*%y + X%*%beta))
	ans <- -1*sum(log(inside))
	return(ans)	
	}

W <- as.matrix(W)

##################################
########## Imputation 1 ##########
##################################

imp1 <- impdata1[1][[1]][[1]]
imp1.models <- list(length=4)

# Model 1
y <- imp1$onset
X <- as.matrix(cbind(1, imp1$warl, imp1$gdpenl, imp1$lpopl1, imp1$lmtnest, imp1$Oil, imp1$nwstate, imp1$instab, imp1$polity2l, imp1$ethfrac, imp1$govpower, imp1$discrim, imp1$powerless, imp1$sepautonomy, imp1$regautonomy, imp1$conc, imp1$relpop, imp1$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M1 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M1.coef <- M1$par
M1.ses <- sqrt(diag(solve(-M1$hessian)))

M1.mat <- cbind(M1.coef, M1.ses)
rownames(M1.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp1.models[[1]] <- M1.mat


# Model 2
y <- imp1$onset
X <- as.matrix(cbind(1, imp1$warl, imp1$gdpenl, imp1$lpopl1, imp1$lmtnest, imp1$Oil, imp1$nwstate, imp1$instab, imp1$anocl, imp1$deml, imp1$ethfrac, imp1$govpower, imp1$discrim, imp1$powerless, imp1$sepautonomy, imp1$regautonomy, imp1$conc, imp1$relpop, imp1$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M2 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M2.coef <- M2$par
M2.ses <- sqrt(diag(solve(-M2$hessian)))

M2.mat <- cbind(M2.coef, M2.ses)
rownames(M2.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp1.models[[2]] <- M2.mat


# Model 3
y <- imp1$onset
X <- as.matrix(cbind(1, imp1$warl, imp1$gdpenl, imp1$lpopl1, imp1$lmtnest, imp1$Oil, imp1$nwstate, imp1$instab, imp1$anocl, imp1$deml, imp1$ethfrac, imp1$govpower, imp1$discrim, imp1$powerless, imp1$sepautonomy, imp1$regautonomy, imp1$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M3 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M3.coef <- M3$par
M3.ses <- sqrt(diag(solve(-M3$hessian)))

M3.mat <- cbind(M3.coef, M3.ses)
rownames(M3.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "peaceyrs", "W")
imp1.models[[3]] <- M3.mat


# Model 4
y <- imp1$onset
X <- as.matrix(cbind(1, imp1$warl, imp1$gdpenl, imp1$lpopl1, imp1$lmtnest, imp1$Oil, imp1$nwstate, imp1$instab, imp1$anocl, imp1$deml, imp1$ethfrac, imp1$conc, imp1$relpop, imp1$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M4 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M4.coef <- M4$par
M4.ses <- sqrt(diag(solve(-M4$hessian)))

M4.mat <- cbind(M4.coef, M4.ses)
rownames(M4.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "conc", "relpop", "peaceyrs", "W")
imp1.models[[4]] <- M4.mat




##################################
########## Imputation 2 ##########
##################################

imp2 <- impdata1[1][[1]][[2]]
imp2.models <- list(length=4)

# Model 1
y <- imp2$onset
X <- as.matrix(cbind(1, imp2$warl, imp2$gdpenl, imp2$lpopl1, imp2$lmtnest, imp2$Oil, imp2$nwstate, imp2$instab, imp2$polity2l, imp2$ethfrac, imp2$govpower, imp2$discrim, imp2$powerless, imp2$sepautonomy, imp2$regautonomy, imp2$conc, imp2$relpop, imp2$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M1 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M1.coef <- M1$par
M1.ses <- sqrt(diag(solve(-M1$hessian)))

M1.mat <- cbind(M1.coef, M1.ses)
rownames(M1.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp2.models[[1]] <- M1.mat


# Model 2
y <- imp2$onset
X <- as.matrix(cbind(1, imp2$warl, imp2$gdpenl, imp2$lpopl1, imp2$lmtnest, imp2$Oil, imp2$nwstate, imp2$instab, imp2$anocl, imp2$deml, imp2$ethfrac, imp2$govpower, imp2$discrim, imp2$powerless, imp2$sepautonomy, imp2$regautonomy, imp2$conc, imp2$relpop, imp2$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M2 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M2.coef <- M2$par
M2.ses <- sqrt(diag(solve(-M2$hessian)))

M2.mat <- cbind(M2.coef, M2.ses)
rownames(M2.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp2.models[[2]] <- M2.mat


# Model 3
y <- imp2$onset
X <- as.matrix(cbind(1, imp2$warl, imp2$gdpenl, imp2$lpopl1, imp2$lmtnest, imp2$Oil, imp2$nwstate, imp2$instab, imp2$anocl, imp2$deml, imp2$ethfrac, imp2$govpower, imp2$discrim, imp2$powerless, imp2$sepautonomy, imp2$regautonomy, imp2$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M3 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M3.coef <- M3$par
M3.ses <- sqrt(diag(solve(-M3$hessian)))

M3.mat <- cbind(M3.coef, M3.ses)
rownames(M3.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "peaceyrs", "W")
imp2.models[[3]] <- M3.mat


# Model 4
y <- imp2$onset
X <- as.matrix(cbind(1, imp2$warl, imp2$gdpenl, imp2$lpopl1, imp2$lmtnest, imp2$Oil, imp2$nwstate, imp2$instab, imp2$anocl, imp2$deml, imp2$ethfrac, imp2$conc, imp2$relpop, imp2$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M4 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M4.coef <- M4$par
M4.ses <- sqrt(diag(solve(-M4$hessian)))

M4.mat <- cbind(M4.coef, M4.ses)
rownames(M4.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "conc", "relpop", "peaceyrs", "W")
imp2.models[[4]] <- M4.mat




##################################
########## Imputation 3 ##########
##################################

imp3 <- impdata1[1][[1]][[3]]
imp3.models <- list(length=4)

# Model 1
y <- imp3$onset
X <- as.matrix(cbind(1, imp3$warl, imp3$gdpenl, imp3$lpopl1, imp3$lmtnest, imp3$Oil, imp3$nwstate, imp3$instab, imp3$polity2l, imp3$ethfrac, imp3$govpower, imp3$discrim, imp3$powerless, imp3$sepautonomy, imp3$regautonomy, imp3$conc, imp3$relpop, imp3$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M1 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M1.coef <- M1$par
M1.ses <- sqrt(diag(solve(-M1$hessian)))

M1.mat <- cbind(M1.coef, M1.ses)
rownames(M1.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp3.models[[1]] <- M1.mat


# Model 2
y <- imp3$onset
X <- as.matrix(cbind(1, imp3$warl, imp3$gdpenl, imp3$lpopl1, imp3$lmtnest, imp3$Oil, imp3$nwstate, imp3$instab, imp3$anocl, imp3$deml, imp3$ethfrac, imp3$govpower, imp3$discrim, imp3$powerless, imp3$sepautonomy, imp3$regautonomy, imp3$conc, imp3$relpop, imp3$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M2 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M2.coef <- M2$par
M2.ses <- sqrt(diag(solve(-M2$hessian)))

M2.mat <- cbind(M2.coef, M2.ses)
rownames(M2.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp3.models[[2]] <- M2.mat


# Model 3
y <- imp3$onset
X <- as.matrix(cbind(1, imp3$warl, imp3$gdpenl, imp3$lpopl1, imp3$lmtnest, imp3$Oil, imp3$nwstate, imp3$instab, imp3$anocl, imp3$deml, imp3$ethfrac, imp3$govpower, imp3$discrim, imp3$powerless, imp3$sepautonomy, imp3$regautonomy, imp3$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M3 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M3.coef <- M3$par
M3.ses <- sqrt(diag(solve(-M3$hessian)))

M3.mat <- cbind(M3.coef, M3.ses)
rownames(M3.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "peaceyrs", "W")
imp3.models[[3]] <- M3.mat


# Model 4
y <- imp3$onset
X <- as.matrix(cbind(1, imp3$warl, imp3$gdpenl, imp3$lpopl1, imp3$lmtnest, imp3$Oil, imp3$nwstate, imp3$instab, imp3$anocl, imp3$deml, imp3$ethfrac, imp3$conc, imp3$relpop, imp3$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M4 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M4.coef <- M4$par
M4.ses <- sqrt(diag(solve(-M4$hessian)))

M4.mat <- cbind(M4.coef, M4.ses)
rownames(M4.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "conc", "relpop", "peaceyrs", "W")
imp3.models[[4]] <- M4.mat




##################################
########## Imputation 4 ##########
##################################

imp4 <- impdata1[1][[1]][[4]]
imp4.models <- list(length=4)

# Model 1
y <- imp4$onset
X <- as.matrix(cbind(1, imp4$warl, imp4$gdpenl, imp4$lpopl1, imp4$lmtnest, imp4$Oil, imp4$nwstate, imp4$instab, imp4$polity2l, imp4$ethfrac, imp4$govpower, imp4$discrim, imp4$powerless, imp4$sepautonomy, imp4$regautonomy, imp4$conc, imp4$relpop, imp4$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M1 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M1.coef <- M1$par
M1.ses <- sqrt(diag(solve(-M1$hessian)))

M1.mat <- cbind(M1.coef, M1.ses)
rownames(M1.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp4.models[[1]] <- M1.mat


# Model 2
y <- imp4$onset
X <- as.matrix(cbind(1, imp4$warl, imp4$gdpenl, imp4$lpopl1, imp4$lmtnest, imp4$Oil, imp4$nwstate, imp4$instab, imp4$anocl, imp4$deml, imp4$ethfrac, imp4$govpower, imp4$discrim, imp4$powerless, imp4$sepautonomy, imp4$regautonomy, imp4$conc, imp4$relpop, imp4$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M2 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M2.coef <- M2$par
M2.ses <- sqrt(diag(solve(-M2$hessian)))

M2.mat <- cbind(M2.coef, M2.ses)
rownames(M2.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp4.models[[2]] <- M2.mat


# Model 3
y <- imp4$onset
X <- as.matrix(cbind(1, imp4$warl, imp4$gdpenl, imp4$lpopl1, imp4$lmtnest, imp4$Oil, imp4$nwstate, imp4$instab, imp4$anocl, imp4$deml, imp4$ethfrac, imp4$govpower, imp4$discrim, imp4$powerless, imp4$sepautonomy, imp4$regautonomy, imp4$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M3 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M3.coef <- M3$par
M3.ses <- sqrt(diag(solve(-M3$hessian)))

M3.mat <- cbind(M3.coef, M3.ses)
rownames(M3.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "peaceyrs", "W")
imp4.models[[3]] <- M3.mat


# Model 4
y <- imp4$onset
X <- as.matrix(cbind(1, imp4$warl, imp4$gdpenl, imp4$lpopl1, imp4$lmtnest, imp4$Oil, imp4$nwstate, imp4$instab, imp4$anocl, imp4$deml, imp4$ethfrac, imp4$conc, imp4$relpop, imp4$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M4 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M4.coef <- M4$par
M4.ses <- sqrt(diag(solve(-M4$hessian)))

M4.mat <- cbind(M4.coef, M4.ses)
rownames(M4.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "conc", "relpop", "peaceyrs", "W")
imp4.models[[4]] <- M4.mat




##################################
########## Imputation 5 ##########
##################################

imp5 <- impdata1[1][[1]][[5]]
imp5.models <- list(length=4)

# Model 1
y <- imp5$onset
X <- as.matrix(cbind(1, imp5$warl, imp5$gdpenl, imp5$lpopl1, imp5$lmtnest, imp5$Oil, imp5$nwstate, imp5$instab, imp5$polity2l, imp5$ethfrac, imp5$govpower, imp5$discrim, imp5$powerless, imp5$sepautonomy, imp5$regautonomy, imp5$conc, imp5$relpop, imp5$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M1 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M1.coef <- M1$par
M1.ses <- sqrt(diag(solve(-M1$hessian)))

M1.mat <- cbind(M1.coef, M1.ses)
rownames(M1.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp5.models[[1]] <- M1.mat


# Model 2
y <- imp5$onset
X <- as.matrix(cbind(1, imp5$warl, imp5$gdpenl, imp5$lpopl1, imp5$lmtnest, imp5$Oil, imp5$nwstate, imp5$instab, imp5$anocl, imp5$deml, imp5$ethfrac, imp5$govpower, imp5$discrim, imp5$powerless, imp5$sepautonomy, imp5$regautonomy, imp5$conc, imp5$relpop, imp5$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M2 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M2.coef <- M2$par
M2.ses <- sqrt(diag(solve(-M2$hessian)))

M2.mat <- cbind(M2.coef, M2.ses)
rownames(M2.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "conc", "relpop", "peaceyrs", "W")
imp5.models[[2]] <- M2.mat


# Model 3
y <- imp5$onset
X <- as.matrix(cbind(1, imp5$warl, imp5$gdpenl, imp5$lpopl1, imp5$lmtnest, imp5$Oil, imp5$nwstate, imp5$instab, imp5$anocl, imp5$deml, imp5$ethfrac, imp5$govpower, imp5$discrim, imp5$powerless, imp5$sepautonomy, imp5$regautonomy, imp5$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M3 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M3.coef <- M3$par
M3.ses <- sqrt(diag(solve(-M3$hessian)))

M3.mat <- cbind(M3.coef, M3.ses)
rownames(M3.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "govpower", "discrim", "powerless", "sepautonomy", "regautonomy", "peaceyrs", "W")
imp5.models[[3]] <- M3.mat


# Model 4
y <- imp5$onset
X <- as.matrix(cbind(1, imp5$warl, imp5$gdpenl, imp5$lpopl1, imp5$lmtnest, imp5$Oil, imp5$nwstate, imp5$instab, imp5$anocl, imp5$deml, imp5$ethfrac, imp5$conc, imp5$relpop, imp5$peaceyrs))
par <- rep(0.1, (ncol(X)+1))

M4 <- optim(par, fn=ll.Wlogit, y=y, W=W, X=X, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
M4.coef <- M4$par
M4.ses <- sqrt(diag(solve(-M4$hessian)))

M4.mat <- cbind(M4.coef, M4.ses)
rownames(M4.mat) <- c("Constant", "warl", "gdpenl", "lpopl1", "lmtnest", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "conc", "relpop", "peaceyrs", "W")
imp5.models[[4]] <- M4.mat

##################################

save(imp1.models, file="imp1M.RData")
save(imp2.models, file="imp2M.RData")
save(imp3.models, file="imp3M.RData")
save(imp4.models, file="imp4M.RData")
save(imp5.models, file="imp5M.RData")

